Absence of spin liquid in non-frustrated correlated systems 
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The question of the existence of a spin hquid state in the half-filled Hubbard model on the hon- 
eycomb (aka graphene) lattice is revisited. The Variational Chister Approximation (VGA), the 
Cluster Dynamical Mean Field Theory (CDMFT) and the Cluster Dynamical Impurity Approxima- 
tion (CDIA) are applied to various cluster systems. Assuming that the spin liquid phase coincides 
with the Mott insulating phase in this non-frustrated system, we find that the Mott transition is 
pre-empted by a magnetic transition occuring at a lower value of the interaction U , and therefore the 
spin liquid phase does not occur. This conclusion is obtained using clusters with two bath orbitals 
connected to each boundary cluster site. We argue that using a single bath orbital per boundary 
site is insufficient and leads to the erroneous conclusion that the system is gapped for all nonzero 
values of U . 



INTRODUCTION 

There is currently a keen interest in materials and mod- 
els that could display a spin liquid state. Such a state is 
characterized by the presence of local magnetic moments 
that do not order at any temperature; spin-spin correla- 
tions then decay exponentially as a function of distance 
(or algebraically, in so-called algebraic spin liquids). In 
theoretical language, it can be described as an insula- 
tor that is not adiabatically connected to a band insula- 
tor, but is rather a pure Mott insulator, without sponta- 
neously broken spatial or spin symmetries. 

A spin liquid could arise from the presence of strong 
geometric frustration, for instance in materials with a 
structure resembling that of the Kagome lattice, such 
as Herbertsmithite, or other kind of frustrated geome- 
tries. It has been conjectured that spin liquids could also 
appear in the intermediate-coupling regime of strongly- 
correlated systems, somewhere between a metallic (or 
semi-metallic) phase and a magnetic phase. It is the 
latter possibility that we will entertain in this work. 

More specifically, we will address the controversy about 
the existence of a spin liquid in the phase diagram of the 
half-filled Hubbard model on the honeycomb lattice. The 
corresponding Hamiltonian reads 

H^-tY, {4aCja + H.c.) +UJ2 "*T"4 (1) 

(ij) i 

where Cia annihilates a fermion of spin projection a =ti i 
at site i, n„ = c^Ca is the number of fermions of spin a 
at site I, and (ij) denotes the nearest-neighbor pairs on 
the lattice. This model attempts minimally at describing 
electron-electron interactions in a material like graphene, 
although a realistic description of graphene should in- 
volve long-range Coulomb interactions, phonons, and so 
on. On the other hand, systems of ultracold atoms 
could be arranged to be fairly accurately described by 
the Hamilonian ([I]), with adjustable interaction strength 
U. 



In a recent work on the matter, Meng et aL[T], us- 
ing Quantum Monte Carlo simulations, have argued that 
a spin liquid phase exists in Model ([T]) in the range 
3.5 < U < 4.3. Below U ~ 3.5, the model is in a semi- 
metallic state, and beyond U ~ 4.3 is is in a antifer- 
romagnetic state (the two sublattices carrying opposite 
magnetizations). But more recently, this result has been 
challenged by Sorella et al. pi , also using Quantum Monte 
Carlo simulations, albeit with larger systems. Although 
this controversy seems a rather technical one, rooted in 
Monte Carlo methods, it also shows that the model in 
question, if not in a spin liquid state, is very close to one. 

The presence of a spin liquid phase has been sup- 
ported by other works [3HS] using quantum cluster meth- 
ods [6j, such as the Variational Cluster Approximation 
( VCA) [71 [8] , the Cluster Dynamical Impurity Approxi- 
mation (CDIA) |H1 [H] and Cluster Dynamical Mean Field 
Theory (CDMFT) [TOHS]- Quantum cluster methods 
have been used extensively in the last decade to refine 
our understanding of the Mott-Hubbard transition and 
of competing orders (magnetism, superconductivity) in 
strongly correlated materials. They are based on some 
representation of the full systems by a small, finite clus- 
ter of sites, with additional uncorrelated orbitals (the 
'bath') and/or adjustable source terms in the Hamilto- 
nian. These additional elements are determined either by 
a self-consistency or by a variational principle. [I4j More 
specifically, Yu et al. [3] have studied the question within 
the Kane-Mele model, which reduces to Model ([l]) in the 
special case A = 0. They support the existence of a 
spin liquid phase in the range 3 < i7 ^ i, based on the 
computation of the spectral function and the associated 
spectral gap within VCA. Wuet al.^ conclude similarly 
with CDMFT using a QVIC solver. Seki and Ohta 
use the VCA and CDIA to study the specific question of 
the antiferromagnetic transition and the metal-insulator 
transition in Model ([T]). They conclude that the single- 
particle gap opens up at an infinitesimal value of U . This 
also supports the existence of a spin liquid state. 
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FIG. 1. (Color online) Clusters used in this paper. The first 
two pertain to the square lattice Hubbard model, the other 
to the honeycomb lattice. Blue squares represent bath sites, 
black circles cluster sites per se. Dashed lines represent inter- 
cluster links when more than one cluster make up the repeated 
unit of the super-lattice. 



In this work we will argue, on the contrary, that 
Model ([I]) does not lead to a spin liquid phase in the 
intermediate coupling regime. Instead, the transition to- 
wards a spin liquid is pre-empted by a magnetic tran- 
sition towards an antiferromagnetic state, like on the 
square lattice. We will also use quantum cluster meth- 
ods with an exact diagonalization solver (mostly CDMFT 
and CDIA), except that larger bath systems will be used. 
Indeed, we assert that probing the Mott transition with 
a bath system of insufficient size may lead to the wrong 
conclusion. 

The square and honeycomb lattice are both bipartite, 
and the half-filled Hubbard model defined on both lat- 
tices benefits from particle-hole symmetry, which sets the 
value of the chemical potential /i to U/2, and imposes 
constraints on the bath parameters used in CDMFT and 
CDIA. We shall therefore start with a discussion of the 
square lattice model, in which the same methodological 
issues arise, in order to put the honeycomb lattice results 
in perspective. 

In the square-lattice Hubbard model at half-filling, it is 
well-known that the Mott transition is pre-empted by the 
onset of antiferromagnetic order. Nevertheless, the Mott 
transition may be investigated by quantum cluster meth- 
ods if AF order is not allowed to set in. This is how it was 
observed in Ref. IHl In that paper, the systems s4-4b and 
s4-8b (see Fig. [T]) were treated with CDIA. Particle-hole 
symmetry left only one independent variational parame- 
ter in s4-4b: the cluster-bath hybridization parameter 0. 
In System s4-8b, a bath energy ±e was also introduced, 
with opposite signs on the two bath orbitals linked to 
the same cluster site, so as to reflect particle-hole sym- 
metry. A clear first-order Mott transition, with metallic, 
insulating, and unstable solutions, was found for the two 
systems, as shown in Fig. [2] The values of Uci and Uc2 
are shown on these figures by vertical dashed lines. 

To ascertain the presence of a spectral gap, we proceed 
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FIG. 2. (Color online) Full line: Hybridization parameter 9 
as a function of U, obtained in CDIA for System s4-4b (top) 
and s4-8b (bottom). The three solutions (metallic, unstable, 
insulating) are shown. The vertical dashed lines indicate Uci 
and Uc2- Dashed curve: density of states A'^(O). Note the 
non-vanishing DoS in the "insulating" solution of s4-4b. 



as follows: The density of states N{uj) can be computed 
by numerically integrating the spectral function A(k,w) 
over wavevectors. We compute N{uj + irj) at a; = for a 
few values of the Lorenzian broadening rj and extrapolate 
?7 — > using a polynomial fit. The result of this extrapo- 
lation should vanish in the insulating solution, but not in 
the metallic solution. This extrapolated density of states 
is shown (dashed curves) in Fig. |2] The remarkable fea- 
ture is that is does not vanish in the insulating solution 
associated with the system s4-4b, but does, as it should, 
in the large bath system s4-8b. Thus, even though Sys- 
tem s4-4b displays the a first-order transition that has all 
the appearances of a Mott transition, its spectral func- 
tion in the strong-coupling phase has no gap, and thus 
this system does not adequately describe a Mott insula- 
tor. In the context of CDMFT or CDIA, this is related to 
the presence of a single bath orbital per cluster site, and 
to particle-hole symmetry. The latter forces the bath en- 
ergy £ to vanish. A correct description of the insulating 
state requires rather a minimum of two bath orbitals per 
cluster site, with equal and opposite bath energies ±e. 

The presence or not of a gap may also be ascertained by 
computing the particle density n(/Lt) around the particle- 
hole symmetric point /U = U/2 and to look for a plateau 
in /i, which would be the signature of a gap (the con- 
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straints on bath parameters stemming from particle-hole 
symmetry must then be released) . This method does not 
require a Lorenzian broadening, as it involves integrals 
carried along the imaginary frequency axis. We have 
used both this method and the extrapolation method de- 
scribed above and find the same conclusions (see supple- 
mentary material). 

The lesson to carry from the square lattice is that a 
sound description of the Mott transition (hence of a pu- 
tative spin liquid state) is to be found in a cluster system 
with minimally two bath orbitals per cluster site. 
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FIG. 3. (Color online) Result of VGA calculations on the 6- 
site cluster h6. The order parameter (red continuous curve) 
vanishes below U — 3.75. Blue dashed curve: spectral gap 
A computed without variational parameters (CPT). Green 
dotted curve: A computed from the VGA solution {t' used as 
variational parameter) . 



Let us now turn to the honeycomb lattice. Note that 
the antiferromagnetic order, on that lattice, is at zero 
wavevector. Thus the nesting condition is always satis- 
fied and we expect the AF state to be fully gapped. This 
is observed in all systems studied (e.g., by inspecting the 
spectral function). 

The first system studied is a 6-site cluster (h6 on 
Fig.[T]). It has been treated with VGA, using the nearest- 
neighbor hopping parameter <' appearing in the cluster 
Hamiltonian H' and a staggered magnetization M as 
variational parameters (see supplementary material for 
a brief summary of the method). As shown on Fig. |3j 
the system develops a nonzero staggered magnetization 
Ai ioi U > Un — 3.75i. This value of Un is remarkably 
close to the one found in Ref. [1(3.75 <Un < 3.8). This 
is most likely a happy coincidence, as Un will depend 
on cluster size. However, the spectral gap obtained from 
the poles of the Green function is nonzero for all values 
of U , even those below Un- This agrees with Ref.iSj This 
seems a signature of a spin liquid state, but, as we will 
see below, systems that are better equiped to describe the 
Mott transition will lead us to the opposite conclusion. 
Note that the gap computed from the VGA solution lies 
below the one obtained from the Green function without 
variational parameters, hinting that a better variational 



solution prefers a smaller gap. 
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FIG. 4. (Color online) Top: Solution of the h2-4b system. 
Bath energy e obtained in CDIA ([/ci and (7c2 are indicated 
by dotted vertical lines) and CDMFT. The staggered magne- 
tization M obtained by CDIA is also shown. Bottom: same, 
for System h4-6b. The value [7^ where the metallic and insu- 
lating solutions have the same energy, is also indicated. 

System h6-6b, with one bath orbital per cluster site, 
was also studied, and our calculations agree with those 
of Ref. ,5j the system has a spectral gap for all values 
of U and displays no Mott-insulator transition (see sup- 
plementary material). However, we assert that probing 
the Mott transition in this system is unreliable, just as 
it is in System s4-4b for the square lattice. It is safer to 
use systems with two bath orbitals per cluster site. The 
simplest such system for the honeycomb lattice is h2-4b 
(Fig.[l]). We studied this system both with GDMFT and 
GDIA. At half-filling, particle-hole symmetry demands 
that the bath energies of the two bath orbitals connected 
to the same cluster site be opposite in value (±e) . In the 
non-magnetic state, the two sites of the cluster (and the 
correponding bath sites) are related by left-right symme- 
try, and therefore only two bath parameters remain: one 
bath energy e and one hybridization parameter 9. The 
GDMFT paramagnetic solution for e is shown on the up- 
per panel of Fig.[4j An upturn in the value of e at around 
U = 5.5 signals the Mott transition. But no hysteresis is 
seen when the interaction U is swept upwards and down- 
wards, which means that GDMFT in this case does not 
detect the first-order character of the Mott transition. 

Things are different when GDIA is applied to the same 
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FIG. 5. (Color online) Top panel: density of states in the 
semi-metallic (SM, red) and insulating (MI, blue) solutions 
obtained from System h4-6b at U — 7, with a Lorenzian 
broadening rj = 0.05t. Bottom panel: density of states 
N(uj — O.Olt + irj), extrapolated to — > for different sys- 
tems, as indicated. Only the semi-metallic solutions of h2-4b 
and h4-6b have a non negligible value. 

system. As shown again on the upper panel of Fig. |4j two 
CDIA solutions are found: a semi-metallic solution when 
U is increased and an insulating solution when U is de- 
creased. Each of these stops, respectively at Uc2 and Ud , 
and coexist in the range [Uci,Uc2\- The first-order char- 
acter of the solution is therefore clearly seen in CDIA. 
Using the same extrapolating method used in the square 
lattice case, it is easily verified that the spectral gap van- 
ishes throughout the semi-metallic solution, whereas it 
is nonzero in the insulating solution (see Fig. [s]) . Note 
that in that case, N{u! + irj) is computed at uj = O.Olt 
instead oi u! = 0, since A^(0) is expected to vanish in the 
semi-metallic phase. 

Also shown on Fig.|4]is the antiferromagnetic order pa- 
rameter obtained if the Weiss field M is added to the list 
of variational parameters; this constitutes in fact a mix- 
ture of CDIA and VCA, since the parameters at play 
within Potthoff's variational approach are both bath- 
related (e and 9) and cluster-related (Af). We find that, 
in this system, the critical interaction strength for the 
onset of magnetic order is C/ ~ 3.4, an even smaller value 
than in the 6-site cluster VCA computation. Thus in this 
system the (continuous) magnetic transition pre-empts 
the Mott transition and no spin liquid occurs. 

The larger cluster system h4-6b (see Fig. [T]) with two 
bath orbitals per site was also studied. In this case, two 4- 
site clusters are necessary to form a repeated unit. Each 



cluster is hybridized to 6 bath orbitals, two on each edge 
site (the central site is not coupled to the bath). Again, 
particle-hole symmetry and rotational symmetry make 
for only two independent bath parameters, e and 9, like 
for the smaller system h2-4b. The solution is shown on 
the lower panel of Fig. |4] When CDMFT is applied, 
the Mott transition appears clearly at J7 « 6.35, but no 
hysteresis is observed. Again, CDIA finds a semi-metallic 
and an insulating (spin liquid) solution, which overlap 
between Ud and Uc2- Their energy f2 are equal at an 
intermediate value Uc (indicated on the figure). Like in 
the case of the system h2-4b, the spectral gap vanishes 
in the semi- metallic phase (Fig. [5]) . If the cluster Weiss 
field M is added to the list of variational parameters, the 
CDIA predicts an antiferromagnetic transition at Upf = 
4.0, which again means that the Mott transition is pre- 
empted. Thus, this larger system also rules out a spin 
liquid phase. 

In conclusion: We have applied various quantum clus- 
ter methods to the Hubbard model on a honeycomb 
lattice, in order to investigate the possible emergence 
of a spin liquid state. We make the hypothesis that 
the spin liquid state that might emerge in a strongly 
correlated system without magnetic frustration, such as 
the one studied here, coincides with the Mott insulating 
state. The Mott transition itself cannot be adequately ac- 
counted for by CPT or VCA: the cluster's environment 
must be described by a bath of uncorrelated orbitals, 
i.e., by a dynamical mean field, and this bath must be 
large enough (two bath orbitals connected to a single 
cluster site). This leaves CDMFT or CDIA as adequate 
cluster methods to study the Mott transition. Two non- 
magnetic solutions (a semi-metal and an insulator, aka 
spin liquid) are found, separated by a first-order transi- 
tion. The CDIA is the better approach, since it reveals 
clearly the three critical values Ud, Uc2 and Uc- The 
spectral functions computed from these solutions show 
the persistence of the Dirac cones up to the Mott tran- 
sition, hence the gapless character of the semi-metallic 
solution. The magnetic solution can also be obtained in 
CDIA, and always appears at a much lower value of U 
than the Mott transition. This leads us to assert that 
a spin liquid (aka Mott insulator) does not exist in this 
system: it is pre-empted by magnetic order. The criti- 
cal value of U at which the magnetic solution appears is 
comparable to what is found in large-scale Monte Carlo 
simulations [U 12]. Since the Mott transition is a rather 
local phenomenon, we argue that increasing the cluster 
size, which we cannot do with our exact diagonalization 
solvers, would not affect the value of Uc to the point of 
changing our conclusion. 15J In fact, increasing the clus- 
ter size would likely shift Uc to a slightly higher value. 
Therefore we believe that our conclusion carries over to 
large clusters. 

We thank G. Baskaran, J. -P. Faye, R. Shankar, P.V. 
Sriluckshmy and A.-M. Tremblay for useful discussions. 
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